
/*********************************************************************
RANDOMIZATION TEST CODE FOR BLORA PROJECT
This program performs one run of a randomization test for the main analysis in the Blora paper.
It is then called by blora_ri.ado and run 10,000 times to obtain randomization inference p-values.
It proceeds as follows:
(1) Captures differences in means (test statistic) for the actual data (seven regressions)
(2) Randomly assigns 1860 observations to four treatment groups blocking on village and sex for each simulation run 
	(see Doc_Blora_Design, Section 5 for details)
(3) Creates treatment assignment variables based on randomized data
(4) Captures betas and standard errors for the randomized data.
**********************************************************************/

capture prog drop blorasim_ri_embed
program define blorasim_ri_embed, rclass
    version 11.1	
		syntax varlist [if]
        gettoken y : varlist
        set more off
		preserve
		
	keep `y' master_id_final kec_name kec_code vil_name vill_code Z* sex* enum block* GZ* taxgrp infogrp 	
	
*************************
*OBTAIN REAL OBSERVED
*************************
  
*Reg 1: Effect of tax treatment in low info environment
  	 qui: reg `y' Z2 if GZ_1v2==1
            return scalar b_y_1r=_b[Z2]
            return scalar se_y_1r=_se[Z2]
            local e(df_r)
            g df_y_1r=`e(df_r)'

*Reg 2: Effect of info treatment in windfall environment
     qui: reg `y' Z3 if GZ_1v3==1
            return scalar b_y_2r=_b[Z3]
            return scalar se_y_2r=_se[Z3]
            local e(df_r)
            g df_y_2r=`e(df_r)'

*Reg 3: Effect of info treatment in tax environment
     qui: reg `y' Z4 if GZ_2v4==1
            return scalar b_y_3r=_b[Z4]
            return scalar se_y_3r=_se[Z4]
            local e(df_r)
            g df_y_3r=`e(df_r)'

*Reg 4: Effect of tax treatment in high info environment
     qui: reg `y' Z4 if GZ_3v4==1
        	return scalar b_y_4r=_b[Z4]
            return scalar se_y_4r=_se[Z4]
            local e(df_r)
            g df_y_4r=`e(df_r)'

*Reg 5: ATE of tax treatment
     qui: reg `y' taxgrp
            return scalar b_y_5r=_b[taxgrp]
            return scalar se_y_5r=_se[taxgrp]
            local e(df_r)
            g df_y_5r=`e(df_r)'

*Reg 6: ATE of info treatment
     qui: reg `y' infogrp
            return scalar b_y_6r=_b[infogrp]
            return scalar se_y_6r=_se[infogrp]
            local e(df_r)
            g df_y_6r=`e(df_r)'
			
*Reg 7: Interaction of tax and info treatments
     qui: xi: reg `y' i.taxgrp*infogrp 
            matrix b=e(b)
            matrix V=e(V)
            scalar b3r=b[1,3]
            scalar varb3r=V[3,3]
        	return scalar b_y_7r=b3r
            return scalar se_y_7r=sqrt(varb3r)
            local e(df_r)
            g df_y_7r=`e(df_r)'

/***********************************************
*RANDOM RE-ASSIGNMENT
Note: The code below reproduces the random assignment while maintaining the strategy for blocking on village and sex 
described in Section 5 of the document blora_design_memo_FINAL.pdf
***********************************************/
	sort block	
	egen blocktag=tag(block) 

	
	set seed 20110110
	
	g rand_block=uniform() if blocktag==1	
	replace rand_block=rand_block[_n-1] if rand_block==.
	g rand_i=uniform()

	egen drop=max(rand_i) if block2=="023141" | block2=="072310" | block2=="087340" , by(block2)
	drop if rand_i==drop

	sort vill_code sex_Z rand_block rand_i
	egen block2_tag=tag(block2)

	bysort block3: g count=_n

	foreach item in 12 13 14 23 24 34 {
		g dyad`item'=uniform() if block2_tag==1 & sex_Z==0
		replace dyad`item'=dyad`item'[_n-1] if block2_tag==0 & sex_Z==0
	}

	egen rmax=rmax(dyad12 dyad13 dyad14 dyad23 dyad24 dyad34) if count<=6 & sex_Z==0
	
	g z=.
		replace z=1 if dyad12==rmax & count==1 & sex_Z==0 | dyad12==rmax & count==3 & sex_Z==0 | dyad12==rmax & count==5 & sex_Z==0
		replace z=2 if dyad12==rmax & count==2 & sex_Z==0 | dyad12==rmax & count==4 & sex_Z==0 | dyad12==rmax & count==6 & sex_Z==0
		replace z=1 if dyad13==rmax & count==1 & sex_Z==0 | dyad13==rmax & count==3 & sex_Z==0 | dyad13==rmax & count==5 & sex_Z==0
		replace z=3 if dyad13==rmax & count==2 & sex_Z==0 | dyad13==rmax & count==4 & sex_Z==0 | dyad13==rmax & count==6 & sex_Z==0
		replace z=1 if dyad14==rmax & count==1 & sex_Z==0 | dyad14==rmax & count==3 & sex_Z==0 | dyad14==rmax & count==5 & sex_Z==0
		replace z=4 if dyad14==rmax & count==2 & sex_Z==0 | dyad14==rmax & count==4 & sex_Z==0 | dyad14==rmax & count==6 & sex_Z==0
		replace z=2 if dyad23==rmax & count==1 & sex_Z==0 | dyad23==rmax & count==3 & sex_Z==0 | dyad23==rmax & count==5 & sex_Z==0
		replace z=3 if dyad23==rmax & count==2 & sex_Z==0 | dyad23==rmax & count==4 & sex_Z==0 | dyad23==rmax & count==6 & sex_Z==0
		replace z=2 if dyad24==rmax & count==1 & sex_Z==0 | dyad24==rmax & count==3 & sex_Z==0 | dyad24==rmax & count==5 & sex_Z==0
		replace z=4 if dyad24==rmax & count==2 & sex_Z==0 | dyad24==rmax & count==4 & sex_Z==0 | dyad24==rmax & count==6 & sex_Z==0
		replace z=3 if dyad34==rmax & count==1 & sex_Z==0 | dyad34==rmax & count==3 & sex_Z==0 | dyad34==rmax & count==5 & sex_Z==0
		replace z=4 if dyad34==rmax & count==2 & sex_Z==0 | dyad34==rmax & count==4 & sex_Z==0 | dyad34==rmax & count==6 & sex_Z==0

	foreach num of numlist 1 2 3 4 {
		g tc`num'=0 if count!=. & sex_Z==0
		replace tc`num'=1 if z==`num' & count!=. & sex_Z==0
		egen sumt`num'=sum(tc`num') if count!=. & sex_Z==0, by(vill_code)
	}

	drop tc*

	egen rmax2= rmax(dyad12 dyad13 dyad14 dyad23 dyad24 dyad34) if count==7 & sumt1!=3 & sumt2!=3 & sumt3!=3 & sumt4!=3 & sex_Z==0| count==8  & sumt1!=3 & sumt2!=3 & sumt3!=3 & sumt4!=3 & sex_Z==0
	egen rmax3= rmax(dyad23 dyad24 dyad34) if sumt1==3 & sumt2<3 & sumt3<3 & sumt4<3 & count==7  & sex_Z==0 | sumt1==3 & sumt2<3 & sumt3<3 & sumt4<3 & count==8 & sex_Z==0 
	egen rmax4= rmax(dyad13 dyad14 dyad34) if sumt2==3 & sumt1<3 & sumt3<3 & sumt4<3 & count==7  & sex_Z==0 | sumt2==3 & sumt1<3 & sumt3<3 & sumt4<3 & count==8 & sex_Z==0 
	egen rmax5= rmax(dyad12 dyad14 dyad24) if sumt3==3 & sumt1<3 & sumt2<3 & sumt4<3 & count==7  & sex_Z==0 | sumt3==3 & sumt1<3 & sumt2<3 & sumt4<3 & count==8 & sex_Z==0 
	egen rmax6= rmax(dyad12 dyad13 dyad23) if sumt4==3 & sumt1<3 & sumt2<3 & sumt3<3 & count==7  & sex_Z==0 | sumt4==3 & sumt1<3 & sumt2<3 & sumt3<3 & count==8 & sex_Z==0 

	egen rmax7= rmax(dyad34) if sumt1==3 & sumt2==3 & count==7 & sex_Z==0 | sumt1==3 & sumt2==3 & count==8 & sex_Z==0
	egen rmax8= rmax(dyad24) if sumt1==3 & sumt3==3 & count==7 & sex_Z==0 | sumt1==3 & sumt3==3 & count==8 & sex_Z==0
	egen rmax9= rmax(dyad23) if sumt1==3 & sumt4==3 & count==7 & sex_Z==0 | sumt1==3 & sumt4==3 & count==8 & sex_Z==0
	egen rmax10=rmax(dyad14) if sumt2==3 & sumt3==3 & count==7 & sex_Z==0 | sumt2==3 & sumt3==3 & count==8 & sex_Z==0
	egen rmax11=rmax(dyad13) if sumt2==3 & sumt4==3 & count==7 & sex_Z==0 | sumt2==3 & sumt4==3 & count==8 & sex_Z==0
	egen rmax12=rmax(dyad12) if sumt3==3 & sumt4==3 & count==7 & sex_Z==0 | sumt3==3 & sumt4==3 & count==8 & sex_Z==0

	foreach num of numlist 2 3 4 5 6 7 8 9 10 11 12 {
		replace rmax=rmax`num' if rmax==. & rmax`num'!=. & sex_Z==0
	}
	
	drop rmax2 rmax3 rmax4 rmax5 rmax6 rmax7 rmax8 rmax9 rmax10 rmax11 rmax12
	
		replace z=1 if dyad12==rmax & count==7 & sex_Z==0 
		replace z=2 if dyad12==rmax & count==8 & sex_Z==0 
		replace z=1 if dyad13==rmax & count==7 & sex_Z==0 
		replace z=3 if dyad13==rmax & count==8 & sex_Z==0 
		replace z=1 if dyad14==rmax & count==7 & sex_Z==0 
		replace z=4 if dyad14==rmax & count==8 & sex_Z==0 
		replace z=2 if dyad23==rmax & count==7 & sex_Z==0 
		replace z=3 if dyad23==rmax & count==8 & sex_Z==0 
		replace z=2 if dyad24==rmax & count==7 & sex_Z==0 
		replace z=4 if dyad24==rmax & count==8 & sex_Z==0 
		replace z=3 if dyad34==rmax & count==7 & sex_Z==0 
		replace z=4 if dyad34==rmax & count==8 & sex_Z==0 

	drop sum*
	
	foreach num of numlist 1 2 3 4 {
		g tc`num'=0 if count!=. & sex_Z==0
		replace tc`num'=1 if z==`num' & count!=. & sex_Z==0
		egen sumt`num'=sum(tc`num') if count!=. & sex_Z==0, by(vill_code)
	}

	drop tc*

	egen rmax2= rmax(dyad12 dyad13 dyad14 dyad23 dyad24 dyad34) if count==9 & sumt1!=3 & sumt2!=3 & sumt3!=3 & sumt4!=3 & sex_Z==0 | count==10  & sumt1!=3 & sumt2!=3 & sumt3!=3 & sumt4!=3 & sex_Z==0
	egen rmax3= rmax(dyad23 dyad24 dyad34) if sumt1==3 & sumt2<3 & sumt3<3 & sumt4<3 & count==9  & sex_Z==0 | sumt1==3 & sumt2<3 & sumt3<3 & sumt4<3 & count==10 & sex_Z==0 
	egen rmax4= rmax(dyad13 dyad14 dyad34) if sumt2==3 & sumt1<3 & sumt3<3 & sumt4<3 & count==9  & sex_Z==0 | sumt2==3 & sumt1<3 & sumt3<3 & sumt4<3 & count==10 & sex_Z==0 
	egen rmax5= rmax(dyad12 dyad14 dyad24) if sumt3==3 & sumt1<3 & sumt2<3 & sumt4<3 & count==9  & sex_Z==0 | sumt3==3 & sumt1<3 & sumt2<3 & sumt4<3 & count==10 & sex_Z==0 
	egen rmax6= rmax(dyad12 dyad13 dyad23) if sumt4==3 & sumt1<3 & sumt2<3 & sumt3<3 & count==9  & sex_Z==0 | sumt4==3 & sumt1<3 & sumt2<3 & sumt3<3 & count==10 & sex_Z==0 
	egen rmax7= rmax(dyad34) if sumt1==3 & sumt2==3 & count==9 & sex_Z==0 | sumt1==3 & sumt2==3 & count==10 & sex_Z==0
	egen rmax8= rmax(dyad24) if sumt1==3 & sumt3==3 & count==9 & sex_Z==0 | sumt1==3 & sumt3==3 & count==10 & sex_Z==0
	egen rmax9= rmax(dyad23) if sumt1==3 & sumt4==3 & count==9 & sex_Z==0 | sumt1==3 & sumt4==3 & count==10 & sex_Z==0
	egen rmax10=rmax(dyad14) if sumt2==3 & sumt3==3 & count==9 & sex_Z==0 | sumt2==3 & sumt3==3 & count==10 & sex_Z==0
	egen rmax11=rmax(dyad13) if sumt2==3 & sumt4==3 & count==9 & sex_Z==0 | sumt2==3 & sumt4==3 & count==10 & sex_Z==0
	egen rmax12=rmax(dyad12) if sumt3==3 & sumt4==3 & count==9 & sex_Z==0 | sumt3==3 & sumt4==3 & count==10 & sex_Z==0

	foreach num of numlist 2 3 4 5 6 7 8 9 10 11 12 {
		replace rmax=rmax`num' if rmax==. & rmax`num'!=. & sex_Z==0
	}
	
	drop rmax2 rmax3 rmax4 rmax5 rmax6 rmax7 rmax8 rmax9 rmax10 rmax11 rmax12
	
	replace z=1 if dyad12==rmax & count==9 & sex_Z==0 
	replace z=2 if dyad12==rmax & count==10 & sex_Z==0 
	replace z=1 if dyad13==rmax & count==9 & sex_Z==0 
	replace z=3 if dyad13==rmax & count==10 & sex_Z==0 
	replace z=1 if dyad14==rmax & count==9 & sex_Z==0 
	replace z=4 if dyad14==rmax & count==10 & sex_Z==0 
	replace z=2 if dyad23==rmax & count==9 & sex_Z==0 
	replace z=3 if dyad23==rmax & count==10 & sex_Z==0 
	replace z=2 if dyad24==rmax & count==9 & sex_Z==0 
	replace z=4 if dyad24==rmax & count==10 & sex_Z==0 
	replace z=3 if dyad34==rmax & count==9 & sex_Z==0 
	replace z=4 if dyad34==rmax & count==10 & sex_Z==0 

*ASSIGN THE MEN
	g dyad=.
		foreach item in 12 13 14 12 24 34 {
			replace dyad=`item' if dyad`item'==rmax & rmax!=.
		}	

	g dyad_men=.
		replace dyad_men=34 if dyad==12
		replace dyad_men=24 if dyad==13
		replace dyad_men=23 if dyad==14
		replace dyad_men=14 if dyad==23
		replace dyad_men=13 if dyad==24
		replace dyad_men=12 if dyad==34

	replace dyad_men=dyad_men[_n-10] if dyad_men==.
	replace dyad=dyad_men if dyad==.
	replace z=floor(dyad/10) if z==. & count==1 & sex_Z==1 | count==3 & sex_Z==1 | count==5 & sex_Z==1 | count==7 & sex_Z==1 | count==9 & sex_Z==1
	replace z=dyad-z[_n-1]*10 if z==. & count==2 & sex_Z==1 | count==4 & sex_Z==1 | count==6 & sex_Z==1 | count==8 & sex_Z==1 | count==10 & sex_Z==1

   	g z1=z==1
    g z2=z==2
    g z3=z==3
    g z4=z==4

    g gz_1v2=z==1 | z==2
    g gz_1v3=z==1 | z==3
    g gz_1v4=z==1 | z==4
    g gz_2v4=z==2 | z==4
    g gz_3v4=z==3 | z==4
    g TAXGRP=z==2 | z==4
    g INFOGRP=z==3| z==4

*************************
*REGRESSIONS ON RE-RANDOMIZED DATA
*************************
     qui: reg `y' z2 if gz_1v2==1
            return scalar b_y_1=_b[z2]
            return scalar se_y_1=_se[z2]
            local e(df_r)
            g df_y_1=`e(df_r)'

     qui: reg `y' z3  if gz_1v3==1
            return scalar b_y_2=_b[z3]
            return scalar se_y_2=_se[z3]
            local e(df_r)
            g df_y_2=`e(df_r)'

     qui: reg `y' z4  if gz_2v4==1
            return scalar b_y_3=_b[z4]
            return scalar se_y_3=_se[z4]
            local e(df_r)
            g df_y_3=`e(df_r)'

     qui: reg `y' z4  if gz_3v4==1
            return scalar b_y_4=_b[z4]
            return scalar se_y_4=_se[z4]
            local e(df_r)
            g df_y_4=`e(df_r)'

     qui: reg `y' TAXGRP 
            return scalar b_y_5=_b[TAXGRP]
            return scalar se_y_5=_se[TAXGRP]
            local e(df_r)
            g df_y_5=`e(df_r)'

     qui: reg `y' INFOGRP 
            return scalar b_y_6=_b[INFOGRP]
            return scalar se_y_6=_se[INFOGRP]
            local e(df_r)
            g df_y_6=`e(df_r)'

     qui: xi: reg `y' i.TAXGRP*INFOGRP 
        	matrix b=e(b)
            matrix V=e(V)
        	scalar b3=b[1,3]
        	scalar varb3=V[3,3]
            return scalar b_y_7=b3
            return scalar se_y_7=sqrt(varb3)
            local e(df_r)
            g df_y_7=`e(df_r)'

restore
end





